Land Cover 2010 to 2020#
From 2010 to 2020, the Denver MSA has maintained control of urban sprawl. While the population increased by 16%, from 2.54 million to 2.95 million (block-level population from the 2010 and 2020 decennial census datasets), the urban area has expanded by only 3.5%, from 1,200 to 1,240 square miles, according to data from MRLC.
Show code cell content
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio as rio
from shapely.geometry import Polygon, Point
crs = "EPSG:2232"
%run assets/colors.py
%run assets/rcParams.py
C:\Users\leeje\OneDrive\G-4\land-use-env-modeling\assignments\sprawl-allocation\memo\assets\rcParams.py:8: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
plt.style.use("seaborn-whitegrid")
We used land cover data from 2008 and 2019 as proxies of the land cover in 2010 and 2020, as these are the best available data from MRLC. Land cover data rasters with 4000*4000 sqft cell sizes. Our study identified 5 types of land cover types:
Urban.
Forests.
Farms.
Wetlands.
Water.
Show code cell content
# Read land cover rasters
# The file paths are relative to the notebook source file's location in the repo's main branch
land_cover_2009 = rio.open("../data/land-cover/2009.tif")
land_cover_2019 = rio.open("../data/land-cover/2019.tif")
def raster_to_shape(raster, crs, to_point=False, value_column="value"):
data_array = raster.read(1)
transform = raster.transform
shapes_data = []
height, width = data_array.shape
if to_point:
for row in range(height):
for col in range(width):
value = data_array[row, col]
if np.isnan(value):
continue
x, y = transform * (col + 0.5, row + 0.5)
# Create a polygon for the cell
point = Point(x, y)
# Add the polygon and its value to the list
shapes_data.append({"geometry": point, value_column: value})
return gpd.GeoDataFrame(shapes_data, crs=crs)
for row in range(height):
for col in range(width):
value = data_array[row, col]
if np.isnan(value):
continue
x_min, y_max = transform * (col, row)
x_max, y_min = transform * (col + 1, row + 1)
# Create a polygon for the cell
polygon = Polygon(
[(x_min, y_min), (x_min, y_max), (x_max, y_max), (x_max, y_min)]
)
# Add the polygon and its value to the list
shapes_data.append({"geometry": polygon, value_column: value})
return gpd.GeoDataFrame(shapes_data, crs=crs)
Show code cell content
# Get a fishnet of info from the land cover rasters
fishnet = (
raster_to_shape(
land_cover_2009, crs, to_point=False, value_column="land_cover_2009"
)
.query("land_cover_2009 != 0")
.copy()
)
points_2019 = (
raster_to_shape(land_cover_2019, crs, to_point=True, value_column="land_cover_2019")
.query("land_cover_2019 != 0")
.copy()
)
# Join points to fishnet, no index right
fishnet = gpd.sjoin(fishnet, points_2019, how="left", predicate="contains").drop(
columns=["index_right"]
)
# Classify land cover types
developed = [21, 22, 23, 24]
land_covers = {
"developed": [21, 22, 23, 24],
"forest": [41, 42, 43],
"farm": [81, 82],
"wetland": [90, 95],
"water": [11],
}
for year in [2009, 2019]:
# First update whether developed bool
fishnet[f"developed_{year}"] = fishnet[f"land_cover_{year}"].isin(developed)
# Then update land cover
fishnet[f"land_cover_type_{year}"] = "other"
for land_cover, codes in land_covers.items():
fishnet.loc[
fishnet[f"land_cover_{year}"].isin(codes), f"land_cover_type_{year}"
] = land_cover
fishnet.drop(columns=[f"land_cover_{year}"], inplace=True)
Show code cell content
import altair as alt
alt.data_transformers.disable_max_rows()
def altair_fishnet(fishnet, column, color_dict, legend_title, title):
chart = alt.Chart(
fishnet.to_crs(4326)[["geometry", column]]
).mark_geoshape().encode(
color=alt.Color(
f"{column}:N",
title=legend_title,
scale=alt.Scale(
domain=list(color_dict.keys()), range=list(color_dict.values())
),
)
).properties(
width=400, height=330, title=title
)
return chart
Show code cell content
from colors import (
palette_primary,
palette_backup2_dimmed,
palette_backup1,
palette_green,
palette_green_faded,
palette_water,
palette_hero,
)
color_dict = {
"developed": palette_hero,
"forest": palette_backup2_dimmed,
"farm": palette_green,
"wetland": palette_green_faded,
"water": palette_water,
"other": "#dddddd",
}
Show code cell source
altair_fishnet(
fishnet=fishnet,
column="land_cover_type_2009",
color_dict=color_dict,
legend_title="Land Cover Type",
title="Land Cover in 2008",
)
Show code cell source
altair_fishnet(
fishnet=fishnet,
column="land_cover_type_2019",
color_dict=color_dict,
legend_title="Land Cover Type",
title="Land Cover in 2019",
)
Show code cell source
fishnet["change"] = "Undeveloped"
fishnet.loc[fishnet["developed_2009"], "change"] = "Developed by 2010"
fishnet.loc[
(fishnet["developed_2019"] & (~fishnet["developed_2009"])), "change"
] = "Developed 2010-2020"
color_dict = {
"Developed by 2010": palette_hero,
"Developed 2010-2020": palette_primary,
"Undeveloped": "#dddddd",
}
altair_fishnet(
fishnet=fishnet,
column="change",
color_dict=color_dict,
legend_title="Development Change",
title="Development Change 2010-2020",
)